Three-dimensional vortex solitons in self-defocusing media 
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We demonstrate that families of vortex solitons are possible in a bi-dispersive three-dimensional 
nonlinear Schrodinger equation. These solutions can be considered as extensions of two-dimensional 
dark vortex solitons which, along the third dimension, remain localized due to the interplay between 
dispersion and nonlinearity. Such vortex solitons can be observed in optical media with normal 
dispersion, normal diffraction, and defocusing nonlinearity 
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Vortex solitons are self-localized solutions of nonlin- 
ear wave equations, which are characterized by a phase 
singularity at the pivotal point. The phase charge of a 
simple closed curve surrounding the vortex core is equal 
to 27rra, where m is the integer vorticity of the solution. 
Vortex solitons have been theoretically predicted in the 
context of superfluids [H, 0]. In those early works, a two- 
dimensional (2D) nonlinear Schrodinger (NLS) equation 
with defocusing nonlinearity was shown to support vortex 
soliton solutions whose intensity vanishes at the vortex 
center and asymptotically approaches a constant value at 
infinity. Such dark vortex solitons were experimentally 
observed in a bulk self-defocusing optical medium [3[. 
The stability of nonlinear vortices depends on the vor- 
ticity number, m. Fundamental vortices with m = 1 are 
energetically favorable (and, as a result stable), which im- 
plies that instabilities of other families of solutions may 
result to the formation of (a set of) fundamental vortices. 
In particular, a 2D dark soliton stripe is unstable to long- 
wave symmetry breaking perturbations, leading to the 
generation of fundamental vortex soliton pairs with op- 
posite vorticities 0, H, 0] • Higher order vortices are also 
unstable and break down into fundamental ones. How- 
ever, in the NLS limit no exponentially growing mode 
exists [7] (the instability may be subexponential) and, 
as a result, multicharged vortices are very long-lived ob- 
jects. Strong instabilities of multicharged vortices can 
be triggered by different mechanisms such as dissipation, 
nonlinearity saturation, or anisotropy [8[. 

Another class of solutions feature ring-shaped inten- 
sity profiles and exist in self- focusing media [9[. Such 
solutions are unstable even in the case of saturable non- 
linearity due to azimuthal instabilities and break down to 
a set of fundamental solitons [10]. Ring vortices can be 
stabilized if the model includes a combination of compet- 
ing (self-focusing and self-defocusing) nonlinear terms, 
such as cubic and quintic [u|, LL2[ or quadratic and self- 
defocusing cubic [13] ones. Azimuthal instabilities can 



also be suppressed by appropriately modulating the am- 
plitude of the solution in the angular direction [14]. Re- 
cently a lot of attention has been attracted to the study 
of vortices in periodic lattices. The lattice can stabilize 
families of ring vortices H, [H, [13, El, Gil by trapping the 
intensity in the lattice potential minima. In a different 
setting, vortex solitons were observed in Bose-Einstein 
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condensates 



In contrast to the plethora of theoretical and experi- 
mental works on vortices in two dimensions, only a few 
works have addressed 3D solutions. Stable toroidal soli- 
tons with vorticity 1 were found in systems with compet- 
ing nonlinearities [H, [22[, and a variety of 3D discrete 
solitons of the vortex type were constructed and explored 



in Refs. [23|, l24j]. In addition, propagation and robust- 



ness of 3D ring optical vortices in the atmosphere was ex- 
amined in Ref. [25[ . Focusing properties of bi-dispersive 
(normal dispersion and normal diffraction) optical sys- 
tems have been studied in [H, [2tJ for self-focusing non- 
linearities. 

The subject of the present paper is to find 3D counter- 
parts of the dark optical vortices which were discovered 
long ago in Refs. [H, 0], and to predict experimental 
conditions necessary for their observation. We demon- 
strate that 3D vortices exist in media with the cubic de- 
focusing nonlinearity and normal GVD (group velocity 
dispersion), i.e., the diffraction-dispersion operator is of 
the hyperbolic type. Optical media realizing this model 

HQ. 



In 



are available, such as specific AlGaAs alloys 
the transverse (x, y) plane, these solutions have the form 
of a dark vortex, whereas along the longitudinal axis 
(i.e., in the temporal direction) they remain localized as 
bright temporal solitons, due to the interplay between 
the normal GVD and defocusing nonlinearity. We first 
construct the solutions in a semi-analytical (and quite 
accurate) form, making use of the Hartree approxima- 
tion [30] - Then, we employ the Newton iteration method 
to find the 3D vortices as numerical solutions to the un- 
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derlying NLS equation. The Hartree approximation pro- 
vides the Newton's method with appropriate initial con- 
ditions. The stability of the vortices is tested by direct 
numerical simulations. We conclude that 3D vortices are 
stable (as long as a super-Gaussian carrying the vortex in 
the transverse plane, which is as a part of the numerical 
procedure, does not suffer essential diffraction). Con- 
sidering applications, the central hole in such a vortex 
pancake may be used as an optically induced aperture to 
dynamically control a probe beam passing through it. 

We start the analysis by introducing the normalized 
NLS equation, 

+ \ (viv - itt) - M 2 ^ = o, (i) 

where Vj_ = <9 2 + <9 2 is the diffraction operator, with 
the normal GVD dispersion and defocusing nonlinearity 
coefficients scaled to unity. In addition to the aforemen- 
tioned optical media (AlGaAs alloys), Eq. (pQ) applies, as 
the Gross-Pitaevskii equation, to Bose-Einstein conden- 
sate in an optical lattice, where z is now time, and t a 
spatial coordinate across the lattice. In this case, the 
negative diffraction along t is achieved when the Bloch 
momentum corresponds to a negative effective mass. 
Notice that if i/j(x, y, z, t) is a solution of Eq. (pQ), then 
= ai(j(ax, ay, at, a 2 z), with real free parameter a, is a 
solution too. Using this scale-invariance property, a one- 
parameter family of solutions can be generated from a 
single solution of Eq. (pQ). Introducing polar coordinates 
(p, 4>) in the (x, y) plane, we look for solutions of Eq. (pQ) 
in the form of 

ip(p, 0, t, z) = u{p, t) ex.-p(—ikz) exp(im^). (2) 

where m is integer vorticity. Substituting Eq. ([2|) in 
Eq. (pQ) and using the above-mentioned scale invariance 
to set k = 1/2, we obtain: 

u + (p~ 1 u p + u pp — m 2 p~ 2 u) — u tt — 2u 3 — 0. (3) 

The asymptotic expansion of the t-independent prob- 
lem at p — > oo yields u = (l/y/2)(l — (m 2 /2)p~ 2 — 
(m 2 /8)(8 + m 2 )p- 4 ) + 0(p~ 6 ). On the other hand, 
the asymptotic expansion at p — > is u = c(pl m l — 
(l/(4(m + l))pl m l+ 2 ) + (9(pl m l+ 4 ). Accordingly, we 
look for solutions of the time-dependent problem with 
lim^o u(p, i) = and lim^^oo u(p, t) = u^i). In the 
latter limit, Eq. (|3]) reduces to the dynamical system as- 
sociated with the single-soliton solution of the NLS equa- 
tion, Uoo — Uoojt — 2^^ = 0, whose commonly known 
soliton solution is 

lim u{p,t) — Uoo(t) — sech(t), (4) 

in compliance with our objective to find solutions that 
look as bright solitons in the temporal direction. 



The total energy conserved by Eq. (pQ) is 

/oo />oo />oo 
/ / \i/j(x,y,t)\ 2 dxdydt. (5) 
-oo J — oo J — oo 

However, this integral of motion diverges for solutions 
with asymptotic form (|4j). A renormalized (convergent) 
form of the energy, that takes into account the asymp- 
totic form of the solution, may be defined as 

/OO POO pOO 
/ / {u1{t)-\^\ 2 )dxdydt. (6) 
-oo J — oo J —oo 

Notice that i/) s = ^oo(t)e^/ 2 = sech(£)e 22; / 2 is a solution 
of Eq. (pQ) representing a bright soliton stripe in three 
dimensions. As a result, both terms in Eq. (|6|) obey the 
above-mentioned scale invariance, which can be used to 
derive a rule for generating a one-parameter family of 
vortex solitons from a single one: P n (a) = P n (l)/a. 

We start the analysis by resorting to the Hartree 
approximation (HA), which is based on the product 
ansatz [30|, 

u(p,t) = R(p)T(t). (7) 

This approximation will then be used as an initial guess 
for numerical solutions based on the Newton's method. 
Substituting Eq. ([7j) in Eq. (|3]), we arrive at a formal 
equation, 

(p- 1 ^ + R" - m 2 p- 2 R + R)T- RT" - R 3 T 3 = (8) 

First, Eq. (|HJ) is to be multiplied with R, integrated 
over p from to pi, and divided by J^ 1 R 2 dp. Taking 
the limit of pi — > oo, the equation 

T — T" — 2T 3 = 0. (9) 

is derived by assuming (without loss of generality) that 
R(p — > oo) = 1. A relevant solution to Eq. © is the 
same temporal bright soliton as the one obtained above 
as the asymptotic wave form, T(t) = sech(t). Next, we 
multiply Eq. (jSj) by T, integrate it from —t\ to £2, divide 
by Jj 2 ^ T 2 (t) dt, and take the limits of ti — > 00, £2 —> 00 ? 
to obtain 

3 (p- 1 ^ + R" - m 2 p- 2 R) + A(R - R 3 ) = 0. (10) 

Equation ([TO]) can be solved numerically by dint of stan- 
dard two-point boundary- value methods, such as shoot- 
ing. 

All our direct numerical simulations are performed us- 
ing the beam propagation method, where the linear part 
is solved using the fast Fourier transform and the non- 
linear part by direct integration (see, for example [3lj). 
Since the solution does not vanish at infinity, and aiming 
to avoid artifacts produced by reflections from domain 
boundaries, we took a 12-th order super-Gaussian (in p) 
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FIG. 1: Shape of the fundamental (m = 1) vortex soliton 
solution, at z = 0, 1, 2, 3. The initial condition is taken as per 
the Hartree approximation. Here and below, contour plots 
display the wave amplitude, versus t and p. 

as a finite-extension carrier for the solution. Of course, 
diffraction of the super-Gaussian cannot be avoided, 
hence, after a large but finite propagation distance, the 
field starts to decay. Eventually, this leads to destabi- 
lization of the vortex. 

In Fig. [1] the evolution of an initial configuration sug- 
gested by the HA is depicted. Notice that the simulation 
window is truncated for large values of p so as to dis- 
play only the evolution of the vortex soliton while elimi- 
nating the above-mentioned irrelevant effect of the back- 
ground diffraction. In all simulations, the initial con- 
ditions did not develop azimuthal instabilities, i.e., the 
intensity distribution corresponding to the vortex soliton 
did not generate any dependence on (/>, nor did it de- 
velop any other instability. Furthermore, as the vortex 
propagates, it slightly broadens in t close to its center. 
This effect does not imply any trend to decay of the vor- 
tex. It is, rather, a manifestation of the relaxation of the 
initial approximate wave form towards an exact vortex 
state. As shown below, this conclusion is in agreement 
with results generated by the Newton iteration method. 

To generate numerically exact stationary vortex- 
soliton solutions, the HA was fed, as an initial guess, 
into the Newton's method for Eq. (j3]). Notice that the 
discretization of Eq. (j3]) in variable p requires special 
attention, since the error is p-dependent. More specif- 
ically, the discretization error becomes larger as p de- 
creases, whereas at large values of p the error is almost 
isotropic (independent of p, to the leading order). The 



FIG. 2: Isosurface plot, corresponding to \u\ 2 = 0.5, of the 
fundamental vortex soliton solution obtained by means of the 
Newton iteration method. 



iteration provided for the convergence of solutions fed by 
the above-mentioned initial HA configurations. The iso- 
intensity profile of the so obtained fundamental vortex 
soliton is depicted in Fig. [2j 

In Fig. El direct simulations of the propagation of a 3D 
vortex soliton, with the initial condition generated by the 
Newton's iteration method, is shown (as above, the ini- 
tial configuration was actually multiplied by the 12-th 
order super-Gaussian). Notice that the vortex soliton 
profile remains invariant in the course of the propaga- 
tion, i.e., the vortex soliton is stable. Comparison with 
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FIG. 3: The same as in Fig. \T\ but if the initial condition 
was taken as a stationary solution generated by the Newton's 
iteration method. 
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z = 



z = 3 




FIG. 4: Propagation of a fundamental (m = 1) vortex soliton 
up to z — 0, 3 at which point the nonlinearity is switched off. 



Fig. [TJ which showed similar evolution initiated by the 
HA, attests to the accuracy of that approximation. A 
noteworthy feature revealed by Fig. [3] (in a more salient 
form than by the HA, cf. Fig.[Tj) is that, close to the vor- 
tex core (at small values of p), the vortex is wider in the 
^-direction. We stress that the vortical phase structure 
of the solution is maintained during the propagation. 

To directly verify that such a 3D vortex soliton solu- 
tion is a nonlinear object indeed, we repeated the same 
simulations, dropping the nonlinear term. As one can see 
in Fig. 2J the vortex quickly diffracts in that case by at 
z = 3. 

We have also performed simulations for double vor- 
tices, with m = 2. They were found to propagate undis- 
torted over large distances, without splitting into funda- 
mental vortices. This observation is, in fact, in agreement 
with a conjecture put forward in Ref. [7[. 

In conclusion, we have demonstrated that vortex soli- 
tons are possible in the bi-dispersive 3D NLS equation. 
These solutions have the form of a dark vortex in the 
spatial plane, whereas they are localized along the tem- 
poral dimension, due to the interplay between the GVD 
and nonlinearity. Such vortex solitons can be observed in 
certain optical media with normal GVD, normal diffrac- 
tion, and defocusing nonlinearity, and in self-repulsive 
BEC. 
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